Lab4

Introduction to spatial data in R

Instructions

In this lab you will begin to get used to using R with spatial data and practice your map making in R for a variety of data types.

Start by reading through sections 7.1-7.3 in Arnold and Tilton.2015. Humanities Data in R.

Then read through ‘Spatial Data Manipulation’ on RSpatial. [http://rspatial.org/spatial/index.html]. There’s a lot here, so feel free to move through it quickly and think of this as reference.

1. From scatterplots to maps: Create a series of point maps with context
2. Map displays: Create a series of simple maps with data
3. Exploring spatial data: Get to know the structure of spatial data in R

Upload to smartsite an Rmarkdown and HTML file.

Due: Monday, April 24, 9a

Results

1. From scatterplots to maps

1. Similar to the steps in section 7.2 of Arnold and Tilton you will construct a map of points on campus, with a context map, using UCplaces.csv and the OpenStreetMap and rgdal libraries.

library (OpenStreetMap)
library (rgdal)

## Loading required package: sp

## rgdal: version: 1.2-5, (SVN revision 648)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-4

datapath = "/users/miadawson/Documents/~UC Davis~/~GEO 200CN/labs/lab 4" #assign directory location to an object called datapath
UCplaces <- read.csv(file.path(datapath, "UCplaces.csv")) #read UC places CSV
map <- openmap (c(38.55, -121.77), c(38.53, -121.74)) #provides the map boundary using upper-left and lower-right corners
mapLatLong <- openproj(map) #converts from the default mercator projection to Latitude Longitude
plot (mapLatLong)
points (UCplaces$long, UCplaces$lat, pch=16) #add the points
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3) #add labels

Create two additional maps: b) change the background context map, [hint: ?openmap will give you a list of possible types]

map2 <- openmap (c(38.55, -121.77), c(38.53, -121.74), type="stamen-toner") #change map type 
mapLatLong2 <- openproj(map2) 
plot (mapLatLong2)
points (UCplaces$long, UCplaces$lat, pch=16) 
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3) 

3. change the extent of the map to include more of Davis, and change the size and appearance of the points and labels as appropriate.

map3 <- openmap (c(38.57, -121.79), c(38.51, -121.71)) #adjust corners to zoom out on map
mapLatLong3 <- openproj(map3)
plot (mapLatLong3)
points (UCplaces$long, UCplaces$lat, pch=5, cex=.5) #change size and shape of points
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3, cex=.8) #change size of labels

2. Map display Using the shapefile ‘State_2010Census_DP1’ which contains 2010 census data for the US, you will follow the steps in 7.3 of Arnold and Tilton to produce:

1. Lat Long map of all the states and Puerto Rico;

library(sp)
library(maptools)

## Checking rgeos availability: TRUE

library(rgeos)

## rgeos version: 0.3-23, (SVN revision 546)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.2-4 
##  Polygon checking: TRUE

state <- readShapeSpatial(fn="State_2010Census_DP1") #read US Census data

## Warning: use rgdal::readOGR or sf::st_read

## Warning: use rgdal::readOGR or sf::st_read

index <- (as.data.frame(state)$STUSPS10 %in% c("AK", "HI")) #identify AK and HI data
state <- state[!index,] #remove AK and HI data
plot(state)

2. a map in the UTM coordinate system with the states labeled [hint: you’ll need to load library(rgeos), and don’t forget to plot(stateTrans), before you use the funtion text() to create your labels.];

#convert projection to UTM
projectionObj <- CRS(projargs="+proj=longlat") 
proj4string(state) <- projectionObj
stateTrans <- spTransform(x=state, CRSobj=CRS("+proj=utm +zone=14"))

#define centroid 
centroid <- gCentroid(spgeom=stateTrans, byid=TRUE)
head(centroid)

## class       : SpatialPoints 
## features    : 6 
## extent      : -1333889, 4104598, 2361966, 5482413  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=14 +ellps=WGS84

#plot map with labels placed at centroid
plot(stateTrans)
text(x=centroid$x, y=centroid$y, label=as.data.frame(stateTrans)$NAME10, cex=0.7)

and c) a cross hatch map showing the proportion of seasonal of recreational housing in each state.

#calculate proportion of housing units (DP0180001) declared as being “seasonal, recreational, or occasional” (DP0180008) 
stateTransData <- as.data.frame(state)
perHouseRec <- stateTransData$DP0180008 / stateTransData$DP0180001

#convert to quantized IDs
bins <- quantile(perHouseRec, seq(0,1,length.out=8))
binId <- findInterval(perHouseRec, bins)
densityVals <- seq_len(length(bins)) * 5
plot(stateTrans, density=densityVals[binId])

You can also create a more familiar choropleth map with solid colors by defining your own color ramp/palette, or using one of the ramps defined within R.

colSet01 = c("grey90", "grey80", "grey70", "grey60", "grey50", "grey40", "grey30", "grey20")
plot (stateTrans, col=colSet01[binId])

colSet02 = rev(heat.colors(length(bins)))
plot (stateTrans, col=colSet02[binId])

4. produce two additional maps using the color palettes defined above as colSet01 and colSet02.

#calculate proportion of population (DP0080001) declared as being “Black or African American alone or in combination with one or more other races ” (DP0090002) 
pctBlack <- stateTransData$DP0090002 / stateTransData$DP0080001

#convert to quantized IDs
bins2 <- quantile(pctBlack, seq(0,1,length.out=8))
binId2 <- findInterval(pctBlack, bins2)
densityVals2 <- seq_len(length(bins2)) * 5
plot (stateTrans, col=colSet01[binId2])

plot (stateTrans, col=colSet02[binId2])

3. Exploring spatial data You will use the shapefile “galap.shp”, which has data on the Galapagos Island including species count, island area and elevation. Using the ‘Spatial Data Manipulation’ on RSpatial. [http://rspatial.org/spatial/index.html] as a reference add code to the Rmarkdown file to complete the following steps and respond to the questions.

For this exercise only use the sp and raster packages (and indirectly, packages such as rgdal that raster may use in the background). Consider answering by writing sentences that use inline R commands like 5+5 = ‘r5+5’. a) Read the shapefile “galap.shp” (in galap.zip) into R to create a SpatialPolygonsDataFrame.

library(raster)
galap <- shapefile("galap.shp") #read galap data

2. Show the class of the object in R.

class(galap)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

The class is a “SpatialPolygonsDataFrame”.

3. How many polygons are there in this object?

length(galap@polygons)

## [1] 30

There are 30 polygons in this object.

4. How many variables are there in this object?

dim(galap) 

## [1] 30 12

There are 12 variables.

5. Make a plot of the islands (their outlines).

plot(galap)

6. Make a scatter plot of the number of species as a function of the size of the island.

plot(galap$species~galap$area)

7. What quantity would you use to make a choropleth to represent the number of species on each island?

I would use species count divided by the log of area. Based on the scatterplot above, I wouldn’t want to simply divide number of species by units of area, because the data is not well distributed. Taking the log of the area normalizes the data, as shown in the scatterplot below.

plot(galap$species~log(galap$area))

8. Create that quantity as a new variable in the SpatialPolyonsDataFrame and plot it with spplot

galap$density = galap$species/log(galap$area)
spplot(galap, 'density')

1. Select the largest island and save it to a new shapefile.

which.max(galap$area) #query which polygon has the largest area

## [1] 2

galap$Island[2] #identify the name of this island

## [1] "Isabela"

Isabela=subset(galap, Island=='Isabela') #subset data to only include this island
plot(Isabela)

writeOGR(Isabela, datapath, 'Isabela_shp', driver='ESRI Shapefile') #write a new shapefile for this island

10. Download elevation data for Ecuador (use the function getData in the raster package).

library(raster)
alt=getData('alt', country='Ecuador') #download altitude data for Ecuador

11. Use the crop function in the raster package and then map the elevation data (add the island outlines).

#change projection of galap to match that of alt
alt@crs #identify coordinate system of alt

## CRS arguments:
##  +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

galap_longlat = spTransform(galap, crs(alt)) #transform galapagos shapefile
galap_crop = crop(alt, galap_longlat) #crop altitute to new galapagos shapefile
plot(galap_crop)
plot(galap_longlat, bg='transparent', add=T) #map elevation data with island outlines

12. Create a contour map of elevation in Ecuador.

contour(alt)